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Abstract 

'We present a finite-difference scheme which solves the Stokes problem in the presence of curvilinear non-conforming 
interfaces and provides second-order accuracy on physical field (velocity, vorticity) and especially on pressure. The gist 
of our method is to rely on the Hclmholtz decomposition of the Stokes equation : the pressure problem is then written 
in an integral form devoid of the spurious sources known to be the cause of numerical boundary layer error in most 
■implementations, leading to a discretization which guarantees a strict enforcement of mass conservation. The ghost 
method is furthermore used to implement the boundary values of pressure and vorticity near curved interfaces. 
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1. Introduction 

Obtaining second-order accuracy on pressure, in nume- 
rical simulations of viscous incompressible flows in the pre- 
sence of rigid interfaces, remains a largely open problem, as 
standard algorithms introduce a numerical boundary layer 
that pollutes the pressure field [1, 2]. These errors, which 
arise in the simplest cases - near flat walls conforming 
•with the underlying mesh in flnite element methods [3] - 
become daunting when interfaces cannot be aligned with 
the computational grid. This latter issue, however, is criti- 
cal to the development of Cartesian grid methods - a quite 
•active field - which seek increased efficiency and flexibility 
by avoiding any remeshing even when dealing with curved 
boundaries. 

\ The origin of pressure errors in direct numerical simu- 
■lations can be examined by simply considering the incom- 
pressible Stokes problem : 



X 
J3 



pvit — — Vp + ^Au + f 
V-u = 

u = Uf, on do, 



(la) 
(lb) 
(Ic) 



on a domain f2, with p and /i the fluid density and vis- 
cosity, u and p, velocity and pressure, f bulk forces, and 
Uf, the velocity boundary condition (BC). (We will use 
throughout the standard indicial notation for partial deri- 
vatives). The key problem is that pressure has no explicit 
expression as a function of others fields, but is determined 
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implicitly within the Stokes problem : it is a Lagrange mul- 
tiplier associated with the condition of incompressibility. 

The strategies developed to cope with this difficulty can 
be sorted into two main families. (I) Either mass conser- 
vation is not strictly enforced - as in fractional step or 
projection methods [4, 5] - but attempt is made to re- 
construct the physical pressure from fields computed at 
intermediate times [3, 6, 7, 8, 10]. (II) Or, the continuity 
condition is strictly enforced via e.g. appropriate polyno- 
mial formulations [11], or by relying on potential [12, 13], 
vorticity [14, 15], or mixed potential- vorticity [19, 20] for- 
mulations. The interest of type (I) methods is that they are 
usually based on the usual "velocity-pressure" formulation 
of the Stokes problem [1, 16], which involves only second 
order derivatives; however, they provide 0(At^/^) [4, 5], 
0{M) [8, 9], or at best 0(At3/2) [3^ ^7, 18] convergence 
for p. Type (II) strategies provide second order accuracy 
on pressure, [12, 19] but rely either on matrix formulations 
of higher rank (the bi-Laplacian of potential formulations) 
or on complex, integral forms of the boundary conditions 
(in pure vorticity methods), which can become quite un- 
tractable with complicated boundary geometries. 

Here, we present a method which both (i) is based on 
the natural velocity and pressure fields, and (ii) guarantees 
a strict enforcement of mass conservation. This is perfor- 
med by constructing the discrete problem - in finite dif- 
ferences - on the basis of a Helmholtz decomposition of 
the Stokes problem (detailed in Sec. 2.1), leading quite 
naturally to a mixed velocity-pressure-potential-vorticity 
formulation. Our solver will be shown to achieve second 
order accuracy on all the physical fields including pressure 
even near curved non-conforming interfaces. 
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As in all numerical schemes of type (I), the critical step 
in our algorithm will be to solve for p. This usually is done 
via the following equation : 

Ap = V-f (2) 

which comes after taking the divergence of (la) and using 
V Au = AV u = 0. This expression, however, brings up 
several problems. 

A first difficulty arises because it is tough to ensure that 
the discrete operators Vx^'^\ and A^'') verify 

the correct commutation relations : in naive implementa- 
tions, V-^'^'A^'^^u is non- vanishing near boundaries, which 
amounts to introducing spurious sources in equation (2). 
This problem is dealt with, on square grids, by using the 
rotational form of the Laplacian [21], but we should antici- 
pate further intricacies when dealing with non-conforming 
boundaries. 

A second, more serious, difficulty arises from the fact 
that, even though pressure verifies the Poisson equation (2), 
it cannot be simply seen as the solution of a plain Pois- 
son (Dirichlet or Neumann) problem [3, 17]. The boun- 
dary condition on pressure indeed has no explicit form. 
The Stokes equation itself does introduce constraints on 
the vector components of Vp on dQ, yet viewing them as 
boundary conditions for the pressure problem raise tre- 
mendous difficulties : 

- setting n • Vp, with n the normal vector to an inter- 
face, is a Neumann condition ; it is associated with 
mass flux through the interface i.e. with incompres- 
sibility 

- setting t • Vp, for any vector t tangent to the in- 
terface, fixes the pressure boundary values up to a 
constant, i.e. sets a Dirichlet condition, and is asso- 
ciated with the no slip boundary condition on velo- 
city. 

The Stokes equation thus sets both a Dirichlet and a Neu- 
mann boundary conditions on the pressure problem. It is a 
property of the Stokes problem that these two conditions 
are compatible, i.e. lead to the same solution for pressure. 
Most pressure-based spatial discretizations, however, in- 
troduce inconsistencies : solving the Poisson equation (2) 
with either boundary condition lead to different solutions 
as originally observed by Gresho et al [22]. Guaranteeing 
that both conditions are simultaneously enforced thus re- 
mains a difficult issue and a recurrent theme in the li- 
terature hinges around the appropriate choice of either 
one [3, 23] : in most cases, the no-slip condition and the 
constraints on mass transport through the interface are 
not simultaneously enforced, which introduces large, un- 
controlled, errors and prevents proper convergence of nu- 
merical approximations for pressure [1]. 

We will show that by performing a Ifclmholtz decom- 
position of the Stokes equation it is possible to formu- 
late the pressure problem in an integral form which gua- 
rantees that its boundary conditions are well-posed and 
consistent with mass conservation especially near bounda- 



ries. To cope with curved boundaries within a finite diffe- 
rence scheme, we will rely on the ghost fluid method [24] , 

which has been developed to take into account situations 
where the solutions to an elliptic problem and their de- 
rivatives have jumps at sharp sub-grid interfaces. In this 
method, interface conditions are implemented through the 
discretized matrix problem itself, in contrast with e.g. im- 
mersed boundary methods where they are represented via 
a set of localized sources. The ghost fluid method was 
first applied to two-phase incompressible fiows [24] , then to 
Poisson problems with Dirichlet or mixed boundary condi- 
tions [25, 26] on fixed rigid interfaces. 

In Section 2 we construct step by step the discrete pres- 
sure problem and analyze its convergence properties. Se- 
veral outstanding questions will then remain regarding the 
incorporation of this solver into the Stokes problem : the 
computation of the boundary conditions for p from the ins- 
tantaneous velocity field ; the formiilation of a consistent 
rotational form of the Laplacian. They are discussed in 
Section 3 and shown to permit the construction of a scheme 
with second-order accTiracy on all fields, inchiding pres- 
sure. Test cases are finally presented in Section 4. 

2. Solving for pressure 

2.1. Helmholtz form of the Stokes problem 

2.1.1. Mass conservation & Helmholtz decomposition 

Our approach is based on the strict enforcement of 
mass conservation, which is achieve by devising the dis- 
cretized problem around the Helmholtz structure of the 
Stokes equation. It is inspired by the effectiveness of using 
the rotational form for the Laplacian [3, 18, 31] : wri- 
ting yuAu = — yttVxo; is innocent from the viewpoint of 
continuum equation, yet it permits to avoid introducing 
spurious sources in the discretization of the Stokes equa- 
tion (la). The reason is that the constraint V-'-'^'' A^'^^u = 
is difficult to enforce numerically (it is non-local, i.e. in- 
volves an extended stencil) , while guaranteeing V- ^''^ V x '•'^^ 
is relatively easy with proper definitions of the curl and 
divergence operators. 

In the same spirit, we seek to write the terms —Vp + f 
in a rotational form, to insure that the computation of the 
pressure field respects incompressibility by construction. 
To do so, let us first recall that in both 2 and 3 dimension, 
the velocity field - being solenoidal - can be written as 
u = VxA, with A a potential vector (we do not fix the 
gauge as yet). Introducing the vorticity field w = Vxu, 
the Stokes equation can therefore be written : 

pVxAt + ij,Vxu> + \/p = f (3) 

which can indeed be seen as a Helmholtz decomposition 
for the field of bulk forces. Equation (3) is the basis of our 
representation of the pressure problem. 
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2.1.2. The 2D case 

As the general treatment of the pressure problem in 
3D is somewhat involved, we focus most of this work on 
the 2D case : the discussion of Section 5 will show that 
the analysis of the Helmholtz form of the Stokes problem 
and the methods developed here directly extrapolate to 
3D problems. 

The 2D velocity field is denoted u = (w, v). As usual, in 
two dimensions, the vorticity uje^ and potential A = tlje^ 
are effectively scalar quantities. The curl operator takes 
two different forms when it is applied on a 2D vector field 
or on a scalar potential : we thus write a; = Vxu = Vx—Uy, 
but u = {ipy, -tpx) = V^tp. 

To better expose the Helmholtz structure of the Stokes 
problem 

pV-^?At = -Vp-/xV^w + f (4) 
we introduce the field 



(5) 



so as to write 

V-L?!> + Vp = f (6) 

This equation reveals some form of "conjugation" between 
p and (j) that is further evidenced by constructing integral 
equations relating these two fields. Let us, for this pur- 
pose, consider a curve F C fi, running from point Ato B; 
with t = {tx,ty) the normalized tangent vector; integra- 
ting along r yields : 



ft 



(7) 



To check that this relation defines p as a univalued func- 
tion, we next need to ascertain that in the case A = i.e. 
when r is a closed loop, the rhs of this equation vanishes. 
Taking T to loop counterclockwise, denoting n = (— <y, tx) 
the outer normal vector, using the Stokes theorem and 
t • W-^cj) = ^ ■ we find : 



- i t-V^(f>+ i f-t = I (-Af^i- Vxf) = 



(8) 



where S in the surface enclosed by F. The last integral 
vanishes because as found after taking the curl of (6) : 
A0 = -Vxf. 

In perfect analogy, path integrals of the form 



j n- (-Vp-Ff) ds 



(9) 



define - up to an irrelevant constant - from p, as in 
particular : 

(10) 



Q = j> n- {-Vp + f) ds 



when F in a closed path. 

These latter integral equations, (9) and (10), are the 
basis of our implementation. By posing the pressure pro- 
blem as a set of path integrals of this form, without ac- 
tually solving for ^, we guarantee that — Vp -h f is the 



curl of an unknown field {4>) ; this in turn guarantees that 
V-(— Vp -|- f) vanishes strictly in the discretized the pro- 
blem. Furthermore, this integral pressure problem is fully 
determined by the boundary values of (p. The solution p of 
this integral problem thus automatically verifies the dual 
Neumann and Dirichlct boundary conditions set by the 
Stokes equation (4) : the agonizing choice [1] of one versus 
another is thus resolved. 

In the rest of this section, we assume that the boundary 
values of (j) are known, and show how (9) and (10) can 
then be discretized using finite differences into a well-posed 
numerical problem for p. 



2.2. The MAC grid 
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Figure 1: The MAC grid 



Our space discretization uses the MAC scheme, which 
involves staggered grids as displayed on Fig. 1, with cell 
size Aa; x Ay. We denote h oc Ax oc Ay a characteristic 
discretization scale. The pressure variables are positioned 
at the center of the cells ( "o" points) ; the first and second 
components of velocity, u and v, on the edges ("A" and 
"o", resp.) ; and the fields ip, u>, and (j) = pi^t + at 
the corners of the cells. These choices guarantee that on 
regular grid points all computations of derivatives only 
require centered differences. 

Integer indices i and j tag the location of the grid lines 
in the x and y directions, respectively : pairs of integer 
indices thus mark the location of the ijj and w-type fields. 
Edges are numbered by the index of their midpoint, of the 
form i + and for horizontal and vertical edges, 

respectively. 

We introduce the differential operator dx, 5y defined by 
their action of discrete fields a"'^ - with here a, (3 integers 
or half-integers : 



{Sxaf^ = 



Ax 
A^ 
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The discrete gradient operator (which normally applies of 
p-type fields) is thus V^'^^ = {5x,6y), the discrete diver- 
gence (which applies on (u, u)-typc fields) has the form 
V-'-'^-'u = dxU + SyV, and the discrete rotational operator 
(which applies of 4>-type fields) is V"""'**^ = {5y,—Sx)- We 
finally denote A^'^) = S'^ + Sy the discrete Laplacian. With 
these definitions, it is clear that 



(11) 



at regular grid points. 



2.3. Paths, edges, and segments 

Our discretization of the pressure problem is construc- 
ted by writing integrals of the form (9) or (10) on all pos- 
sible paths C fl that can be formed using the edges of 
MAC grid cells. These paths can always be decomposed 
as series of rectilinear segments. We define accordingly : 

- a regular segment : any edge of the MAC grid that 
entirely belongs to the fluid domain 

— an irregular segm,ent : the part C of an irregu- 
lar cell edge, i.e. an edge intersected by the domain 
boundary 

Given an edge indexed a, 13, 6°''^ E [0, 1] denotes the frac- 
tion of it that lies within the fluid domain and 7"''^, the 
segment it supports - i.e. its intersection with fl. The seg- 
ments' lengths are thus : 

|y,j±i/2| ^ gi,J±y2^y 

By definition, regular edges {6 = 1) have both ends lying at 
the corner of a cell edge ; irregular edges {0 < 9 < 1) have 
one such regular end point and an irregular one, which lies 
right at the interface ^ . 

Conforming boundaries, which run along the edges of 
MAC grid cells (see Fig. 2), arc treated throughout as a 
special case of non-conforming ones : the edges which be- 
long to the interface are considered as lying outside the 
fluid domain 6 = 0, and the cells they border are, accor- 
dingly, deemed irregular. 

2.3.1. Contour integrals around regular cells 

Let us first consider a regular MAC grid cell such as 
that depicted on Fig. 2-right, and ask how equation (10) for 
the cell contour T = [A, B, C, D, A] can be approximated 
using the discrete values of the fields p and f at neighboring 
MAC grid points. 

The contribution of e.g. the lower edge [A, B] to /p n • 
Vp (with n the outer normal) can be evaluated by Taylor- 
expanding twice (in the x then y direction) around the 



1. We disregard the rare cases when a cell edge would be cros- 
sed several times by the boundary as these become irrelevant for a 
sufficiently small discretization step 
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Figure 2: MAC grid cells near (left) and away from (right) a boun- 
dary 



midpoint M : 
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Ay2 



Ay 



(12) 
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0(Ax^Ay^) 



where e.g. p^j. denotes the yx.x-derivative of the pressure 
field at M. When this expression is added to its coun- 
terpart for the upper edge [C, D] , since —Pyi^ + p'^xx ~ 
0{Ax), the second order terms partly cancel out and their 
difference contributes to third order. Adding the contribu- 
tions from all edges, it finally comes : 



1 



AxAy 



n-Vpds = A'^'^^p+0{Ax^,Ay^) (13) 



As bulk forces are in principle known analytically, we 
could compute directly their contribution to the relevant 
integrals. This however is not too useful in view of the 
errors introduced by the discretization of the pressure gra- 
dients and it is enough to assume, as usual, that the values 
of the X and y components of bulk forces are discretized 
on the 11 and v point (resp.) of the MAC grid. An easy 
treatment for n • f then shows : 

J n • f ds = V-^'^^f + 0(Aa;^ Ay2) (14) 

Comparing equations (13) and (14) leads to an unsur- 
prising expression : 



(15) 



which emphasizes that writing equation (10) around a cell 
contour is a (convoluted) way to derive a discrete approxi- 
mation for the Poisson equation. 

2.3.2. Open path integrals around irregular cells 

The interest of relying on integral expressions unravels 
when the Stokes problem is discretized near irregular MAC 
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grid cells, i.e. cells which intersect the boundary, as depic- 
ted on Fig. 3. We then cannot form closed contours lying 
both on the cell edges and entirely in the fluid domain and 
must hence introduce open paths, i.e. inject in the problem 
certain values of (f) at end points B in the integral ex- 
pressions of the form (9). The only open paths that can 
be used must moreover have their end points lying on the 
interface since this is where the values of cj) are known a 
priori. We will thus use F = [O, B,C, E] in the case of 
Fig. 3-left and [O, C, E] in that of Fig. 3-right. 

We choose not to introduce extra discretization points, 
as is done in other implementations [20] : the pressure field 
is only defined at the centers of MAC grid cells. Estima- 
ting Vp on irregular segments such as [O, B] on Fig. 3-left 
or [O, C] on Fig. 3-right, thus requires using pressure va- 
lues at the centers of the adjoining cells. Some of these 
points may lie beyond the interface : we introduce "ghost" 
pressure points. 



i + I , J (which by definition supports segment 7*+ 2 '-^ ) : 



+ 



pi+ 2 >J+ 2 — 2 '3 2 



(17) 



fyi+2'j 



Ay 



with the same expression modulo x, y-symmetry to define 

€x ^ ■ The above Taylor expansions around the segment 
midpoint showed that : 

et+-^=-(l-^^+-^)^(pf, + /-^) ^^^^ 
+ 0(Ax2,Ay2) 

For the irregular cell depicted on Fig. 3-left, using all 
such expressions to compute the path integral along F = 
[0,B,C,E] now leads to : 





Figure 3; Irregular MAC grid cells 
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(19) 



In the situation illustrated on Fig. 3-right, two Tay- 
lor expansions (in the x and y direction) of p around the 
midpoint M of edge [A, B] lead to the following estimate : 



^j\yix)dx = e'+i'^^^ 



■P' 



Ay 



+ 0(Ax2,Ay2) 



M 

2 ^^v 



(16) 



2.4- The discrete pressure problem 
2.4-1- General form and well-posedness 

In our discretization of the pressure problem, the unk- 
nowns are the values of p on the set composed of the 
centers of all regular and irregular cells. An integral equa- 
tion of the form (9) or (10) is constructed on the contour 
of each of these cells, thus guaranteeing that the number of 
equations equals the number of unknowns. Expression (19) 
can be provided a compact form after defining : 



which reveals that an error term of order 1 arises as soon 
as the integration is performed on an irregular {6 ^ 0, 1) 
cell edge. The integral / n • f can be analyzed along quite 
similar lines. 

2.3.3. Discrete error fields 

To lay out the groundwork for our upcoming conver- 
gence analysis we define, for each (regular or irregular) 
edge of the MAC grid, variables e^'""*^^^^ and e^'^^^^''^ that 
characterize the discretization error. On a horizontal edge 



PE- 90 . , „ 

-- — - — on irregular cells 
AxAy ^ (20) 

on regular ones 







where O and E are the origin and end points of the curve 
r running counterclockwise between the two intersections 
of an irregular cell with the interface. With 9'V^'^^p = 
{95xP,06yp), all integral equations constructed using all 
closed and open paths built using the edges of MAC grid 
cells can finally be written as : 

10] + V-^'^) (^V^'^V) = V-^'') (Of) - V-(^' {Oe) (21) 
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These equations are so far exact, and dropping the last 
term leads to the following approximation : 

M + v-^'') (ev^'^^p^ = v-^'^) (ef) (22) 

which is our discrctizcd pressure problem : the only unk- 
nowns left are now just the values of p on ^lp. This numeri- 
cal scheme, solves (21) with additional sources = V-^"^^ (^e), 
which hence accounts for all the discretization errors. 

The discrete form (22) of the pressure problem presents 
several niceties. First, it is quite simple to implement, as 
it involves nearly usual discrctizcd forms of the divergence 
operator. Second, it guarantees that the discrete pressure 
problem is well-posed : indeed, the underlying consistency 
condition between sources and boundary conditions (the 
discrete form of /^qIi • f = J^V-i), which is required 
to ensure that the discrete problem has a solution, is by 
construction verified on each cell, hence on the full pro- 
blem. Third, its matrix representation is symmetric, thus 
allowing the use of fast algorithms. 

The 0(1) errors present in V-^''-' {9e) near irregular 
cells, however, seem foreboding and the outstanding is- 
sue is now to analyze their consequences for the pressure 
solution. 

2.4-2. Convergence analysis 

Our examination of the numerical errors in the resul- 
ting pressure approximation borrows from the argument 
developed by Jomaa and Macaskill [27] and Gibou and 
Fedkiw [26] to analyze the discretization of the Poisson- 
Dirichlet problem at irregular cells in similar finite-difference 
schemes. They have observed that despite the introduction 
of an error of order in the discretization of the Lapla- 
cian, the solution for the Dirichlet-Poisson problem could 
still be obtained with second order accuracy. The reason 
is that certain discretization errors can be mapped into 
higher order errors in the boundary condition. 

To disentangle, from the discrete field e = {ex,€y), an 
effective error on the boimdary conditions of the pressure 
problem, we now seek to write : 

M + V- (^V('^) g) = V- (^e) (23) 

where g is a discrete field with values on flp (the same 
points as p) and x defined on fi': , the set of all irregiilar 
end points of irregular segments (these points, of course, 
lie on the boundary). 

Let i}p c Up denote the set of all the centers of regular 
cells. The discrete Poisson problem : 



A('^)g = V •('^^ (^e) 
q = 



on Qp 
on Op \ dp 



(24a) 
(24b) 



defines a unique discrete field q with values on Clp. As (24a) 

is written on regular cells only, all 9 values appearing un- 
der the divergence are actually equal to 1, and A^'^^g = 



V.('') (^O'^id.^q^ at these points. Furthermore (see the dis- 
cussion preceding equation (13)), the source term in the 
above Poisson problem, V ■^'^^ (Be) = V (e) = 0{Ax'^) : 
the solution q is therefore also of order 0{Ax'^). 

On regular cells, the field e — V^'^^q is, by construction, 
divergence-free, hence can be written as a discrete curl. In 
other words, there is a scalar field x defined on the set ft^ 
of all corners of regular cells, such that for each regular 
segment 7 : 



+ 



Xb - XA 



(25) 



with I7I the segment length, A and B (with B > A using 
the ordering naturally inherited from the x and y coordi- 
nates) its two end points. Indeed, in perfect analogy with 
the continuum construction of a potential field, x can be 
defined by considering paths F made of regular segments ; 
for F with end points A and B : 



Xb-Xa = ^7 
7cr 



(26) 



where 7 = ±|7| is a signed scalar accounting for the direc- 
tion in which F runs along segment 7. x is well-defined up 
to a constant because equation (24a) guarantees that for 
any closed loop, the rhs vanishes in the above equation. If 
we now consider two points C and £> at a fixed physical 
distance, the difference xc — Xd can be constructed using 
paths containing 0(1/ Ax) segments ; each term in the sum 
in of order 0(|7|) x 0(e) = 0{Ax^), whence x = 0{Ax'^) 
(it is so far defined on the corners of regular cells). 

Requiring additionally that Eq. (25) holds on all ir- 
regular segments now defines x on f2^. On each irregu- 
lar segment 7', the Dirichlet condition (24b) ensures that 

[V'-d)qy = while eT' = 0{Ax) (this was the most pro- 
blematic error) : the constructed values of x on are 
hence again 0{Ax^). 

With these definitions, it is easy to check that Eq. (23) 
holds. Plugging it into equation (21), now yields an equi- 
valent set of exact equations : 

- Xl + V-^'') (oV^^Xp - q)) = V-^'') {0i) (27) 

which shows that the discretization errors introduced when 
solving (22) can be separated into two contributions : (i) 
the usual bulk error q due to the discretization of the 
Laplacian on regular cells ; (ii) a boimdary error which 
perturbs (p at second order. This latter error contributes 
at most second order errors in the underlying Poisson- 
Dirichlet problem defining (6, hence in the solution of the 
pressure problem. All the discretization errors affecting the 
resulting pressure approximation hence arise at second or- 
der. 

2.5. Numerical validation 

To validate our implementation, we will use throughout 
this paper a test Stokes problem introduced in [3]. The 
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time- varying velocity and pressure fields are chosen to be : 

u{x,y) — TT sin(<) sin^(7ra;) sin(27ry) (28a) 
v{x,y) = TT sin(<) sin(27ra::) sin^(7ry) (28b) 
p{x, y) = 47r^ sin(t) cos(7rx) sa\{-ny) (28c) 

and the body forces are computed to match Eq. (1). We 
use the integration domain described in Fig.4-left : 57 = 
[0.4, 2.4]^ \/ where / is an inclusion of radius r — loca- 
ted near the domain center. Using [0.4, 2.4]^ as basis of the 
integration domain permits avoiding trivial error cancella- 
tions due to the symmetry of the solution. The inclusion / 
is also slightly off-center with respect to [0.4, 2.4]^ so as to 
avoid other possible error cancellations due to symmetries 
with respect to the grid. 

The convergence of our solver for the pressure problem 
is tested using (28) at time t = 1. The analytical boundary 
values are then : 

(l){x,y) = 

27r^sin(l) (cos(27ra;)sin^(7ry) -f- cos(27ry)sin^(7rx)) (29) 
— cos(l)sin(27ry)sin^(7ra;) 

Several resolutions x Ny are used, with = Ny = 
20, 40, 80, 160, 320, 640, corresponding to discretization steps 
Ax = Lx/Nx and Ay = Ly/Ny^ with cell dimensions 
Lx — Ly — 2. The matrix representing the discrete pro- 
blem is sparse, allowing us to use band storage ; it is also 
symmetrical, but we have nevertheless used a bi-conjugate 
gradient stabilized algorithm (BICGSTAB) preconditio- 
ned by ILU factorization which we had been using throu- 
ghout this study for its robustness. 

The test results presented in figure 4-right clearly show 
second order for i^, and norms, thus validating 
the convergence analysis performed in the previous section. 









QOL~-norm 




L*JH -norm 




^^l'" -norm 




- ■ slope 2 



-2 -1,5 



Figure 4: Left : The domain used for our numerical test : a square 
domain f! = [0.4,2.4]^ with a cylinder of radius 1/6. Right : the 
results of our convergence analysis. 



time-integration as simple as possible and will consider 
only explicit schemes. Higher-order-in-time and implicit 
schemes based on the ideas developed in the present paper 
are left for future works. 

With Ai the time-step, and denoting it", w", etc. the 
values of hydrodynamic fields at t" = n At, the semi- 
discrete Euler scheme for the Stokes problem simply reads : 



U" 



-vp-MV^w" + f(r) 



At 

V-u"+i = 

u"+i = u(t"+i) on dn 



(30a) 

(30b) 
(30c) 



where we note u(t"+^) (resp. f(i")) the (exact) boundary 
values of the velocity (resp. bulk forces) to be distinguished 
from the approximated solution u"+^. 

To fully exploit our above analysis, we write the cor- 
responding pressure problem as : 



n+l\ 



■0" 



At 



flUJ 



vp - f (r) (31) 



The notation ip{t'^^^) emphasizes that only exactly known 
boundary values are introduced via this term : indeed, 
from the above analysis (see the structure of Eq. (21) 
and (22), or the integral form (9) they derive from) we 
can anticipate that all terms under the operator enter 
the pressure problem only via their boundary values. We 
note that accordingly, boundary values of both ^/)" and 
will have to be estimated from the discrete velocity field 
u". 

3.1. Spatial discretization 

Following Section 2 the pressure problem is discretized 



^ijj{t'^+^) - il)''- 



At 



- jlUJ 



(32) 



while introducing only second order errors. The term |i/)"]/At 
and its relationship with the velocity field can immediately 
be clarified by introducing the discrete currents u", loca- 
lized at the same MAC grid points as the corresponding 
velocities : 



1 



1 



u^^{x^ y) dx 

2 

, z;"(a;,y)dy 



3. Implementation of the Stokes problem 

We now turn to the actual implementation of the uns- 
teady Stokes problem (1) : as our focus remains on the 
details of the spatial discretization, we prefer to keep the 



This definition justifies the introduction of u" values at 
"ghost" grid points, beyond the boundary, where needed 
to account for the mass current passing through irregular 
segments. 

With these definitions, the following relation : 



IVj"1 = v-('" (6*0") 



(33) 
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is exact, and equation (32) is immediately recast as : 



p^(r+i) 



At 



+V ( 



id) 



pu'" 
'At 



(34) 

This discrete problem approximates the set of all integral 
equations (9) around the cells' contours, and requires the 
boundary values of the fields under |.] to be estimated 
at second order on all the intersections of the MAC grid 
with the interface. It provides the values ofp at the centers 
of all regular and irregular cells, which may include a few 
ghost points in the latter case. Accordingly, the gradient of 
pressure as appearing on the rhs of (30a), can be estimated 
from the p values, using only centered differences, at all the 
velocity points that lie on either regular or irregular edges 
(these again, may include a few ghost velocity points). 

Our spatial discretization of the Stokes problem (30) 
relies on the following principles : 

- We only use values of the fields u, v, and p at their 
usual grid points, i.e. at the centers of resp. edges 
and cells. 

- We only use centered differences to estimate the rhs 
of Eq. (30) so as to introduce only second order errors 
in the evaluation of these terms. 

We recognize that introducing extra discretization points 
or more extended stencils in the estimation of these terms 
might be interesting strategies, but it is quite different in 
spirit from what we are developing here, and we therefore 
restrict to the above rules. The second requirement implies 
that to estimate the term W-^uj, we need to define the 
vorticity field on all regular end points of all edges on which 
the velocity field is defined. These again may involve a few 
ghost values. 

Assuming that the values of uj and p are defined on the 
specified points, it is clear that Eq. (30a) can be iterated 
once. With this convention for the discrete points where 

is defined, we can further rewrite equation (34) as : 



At 



id) 



e f 



pu" 
~At 



(35) 



indeed, it is easy to check that evaluating the boundary 
values of w" in (34) by standard (second order) linear 
interpolation on irregular segments, is strictly equivalent 
(numerically) to introducing the term V-^'*-' (BSJ-^lu"). 

Our discretized numerical problem is finally subsumed 
as Eq. (30a) and (35). The remaining difficult questions 
are : given u and v on regular and irregular cell edges at 
time t", how can the corresponding pressure problem be 
defined ? how can the value of uj be determined at all ap- 
propriate points ? how can the discrete fluxes u" be com- 
puted at the required order of approximation (and what is 
the required order of approximation) ? 



o 



a' I 




Figure 5: Localization of tlie velocity and vorticity fields : left, near 
a regular interface ; center, at a regular grid point ; right, near a 
irregular interface. 



3.2. Vorticity estimates 

As usual, the vorticity u — Uy — Vx at a regular point 
can be estimated using centered differences (see illustra- 
tion on Fig. 35-center) : 

oj = Vx'^'^^u + 0{Ax^) = 5tv-5+u + 0{Ax^) (36) 

Let us now consider the case of a regular interface, as 
depicted on Fig. 5-left. In this example of a vertical wall, u 
and all its y-derivatives, especially Uy, are known along the 
boundary : estimating lo at these points only requires to 
introduce an approximation for Vx- This is performed by 
Taylor expanding v in the x direction : taking the abscissa 
of the vertical wall as origin of x coordinates and i indices, 
it comes : 

^ (-8w"J+9w^'^' +0(Aa;2) (37) 



3Ax 



This expression involves only discrete values of v plus its 
boundary value . 

The resulting second order approximations^ for lu = 
Uy — Vx play in our implementation a role analogous to 
Thom's formula and its many variants [1, 28]. However, as 
we do not work in a vorticity-potential formulation, it is 
the values of the velocity field which appear in the rhs. 

We illustrate on Fig. 5-right several situations which 
are encountered when evaluating u near an irregular boun- 
dary : 

- A lies within VL ; 

- B lies on the interface ; 

- C is a ghost point. 

At all these points, uj must be estimated at second order. 
The velocity fields u and v are defined at the centers of 
both regular and irregular segments emanating from A, 
hence oj^ can be simply computed using the regular ex- 
pression (36), which is second order accurate. At point 
B, we can obtain second order accuracy using a classical 
linear interpolation expression involving cj*^ and uj^ , assu- 
ming that uj'^ is second order accurate. The outstanding 
question is how to evaluate uj with second order accuracy 
at all ghost points ? 



2. An alternative way to obtain a second order approximation for 
u) at this regular boundary is to use a linear extrapolation from the 
values of oj : ojO = 2u{d^x) - u){2 /\x) + 0{Ax'^, Ay'^). The sought 
second order accuracy is obtained provided the bulk value of ui are 
known at this order, which can easily be done at regular grid points. 
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Figure 6: Localization of velocity field v and vorticity field to near Using , then M^'^' is immediately obtained by cente- 

an irregular interface on the Mac grid red differences, which finishes to define Ul on ghost points 

such as that illustrated on Fig. 6-right. 



Let us first consider the case depicted on Fig. 6-lcft. 
To avoid singularities due to the possible mismatch bet- 
ween the discrete value i;*+2'J and the analytical boundary 

value need to construct approximations 

which involve the velocity at points that are away from 
the interface. It is convenient to write a general expression 
based on points {i + n+ |, j), and (i + n + |, j). Using a 
Taylor expansion of v along the direction x with origin at 
the ghost point, we find : 

r]Ax \ ) 
+ 0(Aa;2) 

with : 

Q = 8(n + 1) 

^ = 4((l-^'+5'jy- [n 



(38) 




(39) 



77 = Q (1 



)+/3 



1 



■7 



3 



Wc note that in the case n = 0, this expression is identical 
to Eq.(37). Near an irregular interface, it will be used with 
n = 1. 

This approximation for and its counterpart for Uy 
suffice to define w at the ghost point depicted on Fig. 6- 
left. However, few cases also exist when a ghost point lies 
at the end of a single irrcgtilar segment, as illustrated on 
Fig. 6-right. It this latter situation, it is not possible to 
use (38) to estimate Uy. This is performed by first com- 
puting ghost values of u at the points i,j ± ^ via Taylor 
expansion. There again, the mismatch between discrete 
values u'+^'J='=2 and nearby boundary values '•'^5 
may bring in singularities which we avoid by using discrete 
u values which are away from the interface. We thus use : 

rjAa; V ' 7 



0(Aa;2 



3.3. Discrete fluxes 

The discrete fluxes u" appear in the pressure problem 
via the term V-^"^-* (O^) ■ To understand what order of ap- 
proximation is required in the evaluation of these fluxes, let 
us recall that the discretization errors of the pressure pro- 
blem are entirely contained in a term of the form V-*"^^ (6le) : 
we found in Section 2 that second order accuracy on pres- 
sure is achieved even though e may present first order er- 
rors on irregular segments. Consequently, ^ may present 
first order errors on irregular segments (it must of course 
be second order accurate on regular ones). Within our 
time-explicit scheme, as the CFL condition entails that 
At ~ Ax^, u" must be computed with fourth and third 
order accuracy on resp. regular and irregular segments. 

The previous discussion concerning the extrapolation 
of w" to boundary and ghost points relicis on the evaluation 
of both Uy and at the required points. It turns out 
to be convenient to use these intermediary fields in the 
estimation of the discrete fluxes at the appropriate order, 
via expression of the form : 



aAx 



where : 



+ 13 Ax'^ {vi+^^^ - vi^^) + 0{Ax^ 



(42) 



a ■■ 



6 V / 4 
This finishes to define our discrete problem. 
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(43) 



(40) 



3.4. Numerical validation 

We are now in position to implement the numerical 
scheme defined by equations (30) and (35) modulo the 
above technicalities in the computation of w" and u" near 
and at boundaries. To do this, we implement the time- 
dependent test case described in Section 2.5. Panels (a,c,e) 
on Fig. 7-left present the convergence analysis in norms L?, 
and L°°, when averaged over time during numerical 
integration up to time t = 1. These norms are computed 
using either the set of all regular grid cells (dashed lines) 
or using both regular and irregular grid cells (solids lines). 



which may include a few ghost points. A clear second or- 
der convergence is thus found for all the norms and fields 
considered. A few maps of the error field are also provided 
on the right of Fig. 7, panels (b,d,f). 

To validate our algorithm, we also need to check that 
mass conservation 

is properly upheld. This is done by computing the va- 
lues V-^''^(0u") — or regular and irregular cell (see 
Fig. 7-g), which shows faster than second order conver- 
gence. We also show on the same panel that V- (^?u") — 
on regular cells presents only round-off errors, in 
the range of 10^^. 

The above analysis was based on the Euler scheme (30) 
and (35), which can at most provide first order convergence 
in time. To complete our analysis, we thus also consider 
a standard mid-point scheme based on this Euler step. 
Convergence in time is then tested using ^ — 10^^, p = 
1000, with a fixed space discretization = Ny = 40, i.e. 
Ax = Ay = 5 X 10-2 and by varying At = 0.2, 0.4, 0.8, 1.6. 
The results presented in Fig. 8 show again a clear second 
order convergence. 

4. Test cases 

Jf..l. The Faxen problem 

To further test our algorithm, we now consider the case 
depicted on Fig. 9 of a cylinder of radius a moving at a 
constant speed \5 — U ex along the median axis of an in- 
finite horizontal channel of width 2^. We focus on steady 
state, when several approximations exists for the drag co- 
efficient C : Faxen [29] provides an asymptotic expression 
for C in the limit k = a/£ ^ 1 : 



C{k) 



Fxjk) 



47r 



(44) 
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A. 



= -0.9156892732 
= 1.7243844 
-1.730194 
= 2.405644 
= -4.59131 



(45) 



The reverse limit, a/^ — > 1 has been studied by Bungay 
and Brenner [30] : using e = (1 — fc)/fc, they find : 



C(e) = 97rV2£"^/2 ^ 2AB£-'^ + Q^T^/2e-^/^ 
+ (24C + l2D)e-^ + 2ttV2£-^I'^ + ... 



(46) 



where i3, C, and D are integration constants. 

We here use these asymptotic expansions as bench- 
marks tests for our implementation. To avoid difhculties 
due to the progression of the cylinder with respect to the 



Figure 7: From top to bottom : errors on Lt), u, p, and V-u. Left 
convergence analysis. Right : error maps on a 40x40 system. 



gridding, we assume the grid in fixed in the frame of the 
cylinder : the problem is identically mapped onto the case 
of a fixed cylinder, in a channel with moving walls, where 
u = — U. The pressure gradient along the channel - imple- 
mented via homogeneous bulk forces along a; - is computed 
at each timestep to enforce the condition that the average 
fiuid velocity through is also — — U. 

The drag force is in principle given by a path integral 
involving the contour encircling the cylinder (see Fig. 9 for 
details) : 



Ilnds 



(47) 
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(d) Error on divergence 



Nx X Ny. The CFL condition requires the integration time- 
step to be of order 5.10^"^. Our measurements of the drag 
coefficient, presented on Fig. 10, do converge towards the 
relevant asymptotic approximation in both hmits. 
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Figure 10; Comparison between numerical and asymptotic values 
of the drag coefficient C for a cylinder driven in an infinite channel. 
Left : the Faxen case k = a/ 1 <^ 1. Right : the lubrification regime 
(fc^l) 



Figure 8: Convergence analysis on u), u, p, and V u on a 40x40 
system. 




4-2. Flow through a porous medium 

To further ihustrate the capabihties of our algorithm, 
we have implemented a Stokes flow through an array of 
randomly placed cylinder, as illustrated on Fig. 4.2, which 
is a simplified 2D model of a porous medium. Cuts of the 
pressure field along the x direction show clear singularities 
near the obstacles, showing they are very well captured. 



Figure 9: A cylinder of radius a moving at a constant speed U = 
U along the median axis of an infinite horizontal channel of width 

2e 



with the hydrodynamic stress tensor : 



n 



du dv 
^ dy dx 



dy dx ) 



(48) 



As we focus on steady state - which we reach after a reaso- 
nably short transient - momentum conservation V-II = 
guarantees that F can be identically computed on any 
contour F encircling the cylinder : it enables us to use 
contours based on the regular segments of the grid so as 
to avoid introducing extraneous errors due to the estima- 
tion of n along the boundary and only catch errors coming 
from the simulation proper. We have checked that using 
different contours provides identical results modulo round- 
off errors. 

The simulation is implemented, following Ben Richou e^ 
al [13], using a constant radius a — 0.8mm, and a fixed 
spatial step Aa; — Ay — 2.0 x 10"^, while varying the 
channel width and accordingly the number of grid points 





Figure 11: Flow through a porous media. Left : arrow map of the 
velocity field. Right : pressure profile along a horizontal cut through 
the middle of the simulation cell. 



5. Extension to 3D 

Let us now come back to the case of a 3D flow, as 
discussed at the beginning of Section 2.1.1 : writing the 
velocity field as u = VxA, with A a potential vector, the 
Stokes equation reads : 



p\/xAt + fiS/xuj + Vp — f 



(49) 



which can be viewed as the following Helmholtz decompo- 
sition : 

VxB + Vp = f (50) 
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Figure 12: Scheme 



with : 

B = pAt + tiuj . (51) 

The question now is : what form of boundary condi- 
tions can be used to compute p as the scalar quantity 
entering the Helmhohz decomposition (3) ? To answer it, 
let us first examine how this pressure problem can be gi- 
ven an integral representation. We thus consider a volume 
V G fl, with boundary S, that does not intersect the do- 
main boundary dfl. Denoting n the normal vector to any 
surface element dS, we have : 

/ n • Vp dS* = / V-f dV^ (52) 
Js Jv 

This of course is an integral form of the Poisson problem 
Ap — V-f. Second, we consider a volume V that does in- 
tersect the domain boundary. As sketched on Fig. 12, the 
boundary of V can be decomposed into two compo- 
nents : a surface S' which lies in the interior of and a 
surface S" G dfl. We then find : 

/ n-Vpd5= / V-fdF-/ n- (f - VxB) d^ (53) 
Js' Jv Js" 

These integrals representation corresponds to respecti- 
vely equations (10) and (9) in the 2D case : the problem 
is that here, they involve the boundary values of a vector 
field B - versus a scalar field in 2D - which depends on the 
gauge prescription : we a priori do not know how to com- 
pute B and need to guarantee that any gauge fixing will 
lead to an identical discretization for the pressure problem. 

We first write : 

J n VxBd5== J B dl 

with r" the curve enclosing S", and note that the integral 
on the rhs involves only the tangential components of B. 
To study the influence of gauge flxing, we consider one Ag 
such that u — VxA^ : all valid vector potentials for u 
are of the form A = A^ + V(p. We then use the Coulomb 
gauge, V- A — 0, which amounts to requiring that ip verifies 



Aip — — V- As. This Poisson equation can now be provided 
boundary conditions. The choice of any Dirichlet condition 
of (/? - up to an irrelevant constant ~ amounts to (i) fixing 
the components of Vi^ tangent to the interface plus (ii) the 
additional constraint that for any closed curve F S : 
/p Vv?-dl = 0, with dl the normalized line element along F. 
The second condition guarantees that we can define ip via 
integral equations of the form ips — ^Pa = Vi^ • dl where 
integrals are taken along curves e dfl running from A to 
B. The Dirichlet condition on ip is thus equivalent to the 
prescription of the tangential components of A, provided 
for any closed contour F e dfl they verify : 

J A • dl = J As • dl = ju ndS (54) 

In short, any choice of the tangential components of A can 
be made, provided they are consistent with equation (54), 
i.e. they enforce the correct mass fiuxes through the in- 
terface. In practice, near rigid interfaces, it will always 
be possible to compute analytically such components of a 
valid vector potential, without the need to compute A en- 
tirely. To fix B • dl on the domain boundary, it remains to 
evaluate the tangential components of the vorticity there, 
which like in the 2D case can be performed using deriva- 
tives of the components of the velocity field. 

We thus find that the treatment of 3D fiows should 
follow exactly along similar lines to those presented here 
in 2D. 

6. Conclusion 

We have here shown that by viewing the Stokes equa- 
tion as a Helmholtz decomposition of the field of bulk 
forces, it was possible to construct a discretization based 
on the physical p and u fields, which accurately enforces 
mass conservation. The resulting evaluation of the pres- 
sure field in then devoid of numerical errors due to spu- 
rious mass sources and sinks and has been here shown to 
converge with second order even near boundaries. Our im- 
plementation relies on the introduction of ghost points, yet 
seems to correctly capture singularities of the pressure e.g. 
near obstacles such as in a Darcy flow. 

The present work does not pretend to be complete, 
some limitations can be immediately identified, which may 
be more or less challenging : 

1. we have here focussed on the Stokes fiow : this is 
not a strong problem as implementing the non-linear 
terms of the Navier-Stokes equation does not intro- 
duce significant new technicalities ; 

2. we have relied on time-explicit discretizations so as 
to focus on the difficulties arising from the evaluation 
of the pressure field : extending our method to time- 
implicit scheme might require significant revisions in 
the treatment of the vorticity field, however, and this 
we must leave for future works ; 
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3. finally, we have focussed our discussion on the case of 
2D problems, but have also argued the extension to 
the 3D case should be rather immediate and mostly 

technical. 

We hope, however, that our work can open new routes to 
considering the difficult question of implementing direct 
numerical simulations of fluid flows in the presence of in- 
terfaces. 
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